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ABSTRACT 

We use kinetic simulations of diffusive shock acceleration (DSA) to study 
the time-dependent evolution of plane, quasi-parallel, cosmic-ray (CR) modified 
shocks. Thermal leakage injection of low energy CRs and finite Alfven wave prop- 
agation and dissipation are included. Bohm diffusion as well as the diffusion with 
the power-law momentum dependence are modeled. As long as the acceleration 
time scale to relativistic energies is much shorter than the dynamical evolution 
time scale of the shocks, the precursor and subshock transition approach the 
time-asymptotic state, which depends on the shock sonic and Alfvenic Mach 
numbers and the CR injection efficiency. For the diffusion models we employ, 
the shock precursor structure evolves in an approximately self-similar fashion, 
depending only on the similarity variable, xj (u s t). During this self-similar stage, 
the CR distribution at the subshock maintains a characteristic form as it evolves: 
the sum of two power-laws with the slopes determined by the subshock and total 
compression ratios with an exponential cutoff at the highest accelerated momen- 
tum, p ma x(t). Based on the results of the DSA simulations spanning a range of 
Mach numbers, we suggest functional forms for the shock structure parameters, 
from which the aforementioned form of CR spectrum can be constructed. These 
analytic forms may represent approximate solutions to the DSA problem for as- 
trophysical shocks during the self-similar evolutionary stage as well as during the 
steady-state stage if p max is fixed. 
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Introduction 



Diffusive shock acceleration (DSA) is widely accepted as the primary mechanis m through 



which cosmic rays (CRs) are produced in a variety of astrophysical environments (lBelllll978 



Drurylll983l ; iBlandford &: Eichlerl 119871 ) . The most attractive feature of the DSA theory is 
the simple prediction of the power-law momentum distribution of CRs, f{p) oc p _3o "/(°" _1 ) 
(where a is the shock compression ratio) in the test particle limit. For strong, adiabatic 
gas shocks, this gives a power-law index of 4, which is reasonably close to the observed, 
'universal' index of the CR spectra in many environments. 

However, it was recognized early on, through both analytical and numerical calculations, 
that the DSA can be very efficient and that there are highly nonlinear ba ck-reactions from 
CRs to the underlying flows that modify the spectral form, as well (e.g., iMalkov fc Drury 



200ll . for a review). In such CR modified shocks, the pressure from CRs diffusing upstream 
compresses and decelerates the gas smoothly before it enters the dissipative subshock, creat- 
ing a shock precursor and governing the evolution of the flow velocity in the precursor. On 
the other hand, it is primarily the flow velocity through the precursor and the subshock that 
controls the thermal leakage injection and the DSA of CRs. Hence the dynamical structure 
of the flow and the energy spectrum of CRs must evolve together, influencing each other in 
a self-consistent way. 

It is formation of the precursor that causes the momentum distribution of CRs to deviate 
from the simple test-particle power-law distribution. With a realistic momentum-dependent 
diffusion, n(p), the particles of different momenta, p, experience different compressions, de- 
pending on their diffusion length, I dip) = K(p)/u a (where u s is the shock speed). The particles 
just above the injection momentum (pmj) sample mostly the compression across the subshock 
i<J s ), while those near the highest momentum (p max ) experience the greater, total compres- 
sion across the entire shock structure (at). This leads to the particle distribution function 
that behaves as f (p) oc p-S'WK- 1 ) f or p ^ but flattens gradually to f(p) oc p- 3a t/(o- t -i) 
toward p ~ p max ( iDuffy et al.l I1994T ) . 



Analytic solutions for fip) at the shock have been found in steady-state limits under 
speci al conditions; for example, the case of a constant diffusion coefficient (IDrury et al. 
19821 ) and the ca se of steady-state shocks with a fixed ab ove which particles escape 



from the system flMalkovill997l . Il999l ; lAmato fc Blasil 120051 . 120061 ) . In these treatments, the 
self-consistent solutions involve rather complicated transformations and integral equations, 
so are d i fficult to use in general, although they do provide important insights. In particular, 
Malkovl (119991 ) showed that in highly modified, strong, steady shocks (<7 t ^> 1) with a fixed 
p max , the spectrum of CRs flattens to fip) oc p~ 3 - 5 for nip) oc p a with a > 1/2. He also 
argued that the form of the CR spectrum is universal under these conditions, independent 
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of k(p) and at. In an effort to provide more practical description iBerezhko fc Ellison! (119991 ) 
presented a simple approximate model of the CR spectrum at strong, steady shocks in plane- 
parallel geometry. They adopted a three-element, piece- wise power-law form to represent 
the spectrum at non-relativistic, intermediate, and highly relativistic energies. And they 
demonstrated that this model approximately represents the results of their Monte Carlo 
simulations. 



In iKang fc Jones! (12007! ) (Paper I), from kinetic equation simulations of DSA in plane- 
parallel shocks with the Bohm-like diffusion (k oc p), we showed that the CR injection rate 
and the postshock states approach time-asymptotic values, even as the highest momentum 
Pmax{t) continues to increase with time, and that such shocks then evolve in a "self-similar" 
fashion. We then argued that the nonlinear evolution of the shock structure and the CR 
distribution function in this stage may be described approximately in terms of the similarity 
variables, £ = x/(u s t) and Z = ln(p/p- m j)/ ln[p max (t)/p in j]. Based on the self-similar evo- 
lution, we were able to predict the time-asymptotic value of the CR acceleration efficiency 
as a function of shock Mach number for the assumed models of the thermal leakage injec- 
tion and the wave transportation. In those simulations we assumed that the self-generated 
waves provide scatterings sufficient enough to guarantee the Bohm-like diffusion, and that 
the particles do not escape through either an upper momentum boundary or a free-escape 
spatial boundary. So the CR spectrum extended to ever higher momenta, but at the same 
time the particles with the highest momentum spread over the increasing diffusion length 
scale as Z max cx n(p max )/u s oc p max oc t. We note that in Paper I we considered plane-parallel 
shocks with shock Mach number, 2 < M < 80, propagating into the upstream gas with 
either To = 10 4 K or 10 6 K, since we were interested mainly in cosmic structure formation 
shocks. 

The simplicity of the results in Paper I suggested that it might be possible to obtain an 
approximate analytic expression for the CR spectrum in such shocks, but the simulations 
presented in that paper were not sufficient to address that question. Thus we further carried 
out an extensive set of simulations to explore fully the time-dependent behavior of the CR 
distribution in CR modified shocks with shock Mach numbers M > 10. In this paper, from 
the results of these simulations, we suggest practical analytic expressions that can describe 
the shock structure and the energy spectrum of accelerated particles at evolving CR modified 
shocks in plane-parallel geometry, in which the Bohm-like diffusion is valid. 

In realistic shocks, however, once the diffusion length / max becomes comparable to the 
curvature of shocks, or when the growth of waves generated by the CR streaming instability 
is inefficient, the highest energy particles start to escape from the system before they are 
scattered and advected back through the subshock. In such cases, p max is fixed, and the 
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CR spectrum and the shock structure evolve into steady states. So, for comparison, we 
carried out additional simulations for analogous shocks in which the particles are allowed to 
escape from the system once they are accelerated above an upper momentum boundary, p uh . 
Those shocks achieve true steady states and the shock structure and the CR distribution 
become stationary with forms similar to those maintained during the self-similar stage of 
shock evolution. In this sense, our solution is consistent with the analytic solutions for steady 
state shocks obtained in the previous papers mentioned above. 

In the next section we describe the numerical simulations and results. The approximate 
formula for the CR spectrum will be presented and discussed in §3, followed by a summary 
in §4. We also include an appendix that presents simple analytic and empirical expressions 
that can be used to characterize the dynamical properties of CR modified shocks. 



Numerical Calculations 



2.1. Basic equations 



In our kinetic simulations of DSA, we solve the standard gasdynamic equations with 
the CR pr essure terms in the conservative, Euler i an for m for one-dimensional plane-parallel 
geometry fcang et al.lbooi iKang fc Joneslliooi . l2007t ). 



dp d(up) 
dt dx 



0. 



d(pu) | d(pu 2 + P g + P c 
dt dx 



0, 



d(pe g ) | d 
dt 



(pCgU + PgU) 



dPr 



-u- 



+ W(x,t)-L(x,t), 



(2) 
(3) 



dx dx 
where P g and P c are the gas and CR pressures, respectively, e g = Pgl[p{^ g — 1)] + u 2 /2 is 
the total gas energy per unit mass. The remaining variables, except for L and W, have 
the usual meanings. The injection energy loss term, L(x,t), accounts for the energy carried 
away by the suprathermal particles injected into the CR component at the subshock and is 
subtracted from the postshock gas immediately behind the subshock. The gas heating due 
to the Alfven wave dissipation in the upstream region is represented by the term 

BP 

W(x,t) = -v A -^, (4) 



dx 



where va = B/y/Anp is the local Alfven speed (Paper I). These equations can be used to 
describe parallel shocks, where the large-scale magnetic field is aligned with the shock normal 
and the pressure contribution from the turbulent magnetic fields can be neglected. 
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The CR population is evolved by solving the diffusion-convection equation for the pitch- 
angle-averaged distribution function, f(x,p,t), in the form, 

ld_ 

3 dx 



do 

— + 

dt 



u + u w )- 



— -4 
dy 



d_ 

dx 



(5) 



where g = p A f and y = \n(p) (]Skillinglll975al ). Here, n(x,p) is the spatial diffusion coefficient. 
The CR population is isotropized with respect to the local Alfvenic wave turbulence, which 
would in general move at a speed u w with respect to the plasma. Since the Alfven waves 
upstream of the subshock are expected to be established by the streaming instability, the 
wave speed is set there to be u w = va- Downstream, it is likely that the Alfvenic turbulence 
is nearly isotropic, so we use u w = there. 

We consider two models for CR diffusion: Bohm diffusion and power-law diffusion, 

„2 



k b 



K 



K 



Po 

p 



(6) 



with a = 0.5 — 1. Hereafter, the momentum is expressed in units of m p c, where m p is the 
proton mass and c is the speed of light. So, k* is a constant of dimensions of length squared 
over time. As in our previous studies, we consider diffusion both without and with a density 
dependence , p n I p\ that is, either v — or v — 1. The latter case quenches the CR acoustic 
instability (IDruryl 1 19841 ) and approximately accounts for the compressive amplification of 
Alfven waves. Since we do not follow explicitly the amplification of magnetic fields due to 
streaming CRs, we simply assume that the field strength scales with compression and so the 
diffusion coefficient scales inversely with density. Bohm-like diffusion is an idealization of 
what is expected in a dynamically evolving CR modified shock. As discussed in §2.3 the dif- 
fusion coefficient, which results from resonant scattering with Alfven waves, varies inversely 
with the intensity of the resonant waves. The wave intensity is expected to be amplified 
as the shock evolves from upstream, ambient values via the streaming instability. Bohm 
diffusion represents the simplest nonlinear limited model for that process. The very highest 
momentum CRs will encounter ambient wave intensities, so perhaps below levels implied by 
Bohm diffusion. The model as sumes that the streaming instabili ty quickly amplifies those 



waves to nonlinear levels (e.g., ISkillingl Il975bl : iLucek fc Belli |2000| ). We label the quantities 



upstream of the shock precursor by the subscript '0', those immediately upstream of the gas 
subshock by '1', and those downstream by '2'. So, p , for example, stands for the density of 
the upstream gas. 

Equations (1), (2), (3), and (5) are simultaneously integrated by the CRASH (Cosmic- 
Ray Acceleration SHock) code. The detailed description of the CRASH code can be found 
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m lKang et al.l (120021 ) and Paper I. Three features of CRASH are important to our discussion 
below. First, CRASH applies an adaptive mesh refinement technique around the subshock. 
So the precursor structure is adequately resolved to couple the gas to the CRs of low mo- 
menta, whose diffusion lengths can be at least several orders of magnitude smaller that 
the precursor width. Second, CRASH uses a subgrid shock tracking; that is, the subshock 
position is followed accurately within a single cell on the finest mesh refinement level. Conse- 
quently, the effective numerical subshock thickness needed to compute the spatial derivatives 
in equation (jS]) is always less than the single cell size of the finest grid. Third, we calculate 
the exact subshock speed at each time step to adjust the rest frame of the simulation, so 
that the subshock is kept inside the same grid cell throughout. These three features enable 
us to obtain good numerical convergence in our solutions with a minimum of computational 
efforts. As shown in Paper I, the CRASH code can obtain reasonably converged dynamical 
solutions even when the grid spacing in the finest refined level is greater than the diffusion 
length of the lowest energy particles {i.e., Ax 8 > ld(j>m]))- This feature allows us to follow 
the particle acceleration for a large dynamic rage of p m ax/Pmj, typically, ~ 10 9 , although the 
evolution of the energy spectrum at low energies and the early dynamical evolution of the 
shock structure may not be calculated accurately. 



2.2. Simulation Set-up 

The injection and acceleration of CRs at shocks depend in general upon various shock 
parameters such as the Mach number, the magnetic field strength and obliquity angle, and 
the strength of the Alfven turbulence responsible for scattering. In this study we focus 
on the relatively simple case of CR proton acceleration at quasi-parallel shocks, which is 
appropriately described by equations (PQ) - (JHJ) . The details of simulation set-up can be found 
in Paper I, and only a few essential features are briefly summarized here. Except for diffusion 
details, the set-up described here is identical to those reported in Paper I. 

As in Paper I, a shock is specified by the upstream gas temperature T and the initial 
Mach number Mq. Two values of To, 10 4 K and 10 6 K, are considered, representing the warm 
photoionized gas and the hot shock-heated gas often found in astrophysical environments, 
respectively. Then the initial shock speed is given as 

u Stl = c Si0 M = 15 km s-M^J M , (7) 

where c S) o is the sound speed of the upstream gas. All the simulations reported in this 
paper have Mq = 10, which is large enough to produce significant CR modification. In 
Paper I we considered a wide range of shock Mach numbers and examined the Mach-number 
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dependence of the evolution of CR modified shocks. The CR injection and acceleration 
efficiencies are determined mainly by the sonic Mach number an d the relative Alfven Mach 



number for shocks with M > 10 flKang et al.ll2002l ; iKangl 120031 ) . On the other hand, they 



depend sensitively on other model parameters for shocks with lower Mach numbers. In this 
paper we thus focus on the evolution of the CR spectrum at moderately strong shocks with 
^ 10. We will consider the more complicated problem of weaker shocks in a separate 
paper. 

In our problem, three normalization units are required for length, time, and mass. While 
ordinary, one-dimensional, ideal gasdynamic problems do not contain any intrinsic scales, the 
diffusion in the DSA problem introduces one; that is, either a diffusion length or a diffusion 
time, which of course depend on the particle momentum. So let p^ be a specific value of 
the highest momentum that we aim to achieve by the termination time of our simulations. 
Then the greatest width of the precursor is set by the diffusion length of the particles with 
p\ Idip^) = K (.Po,P^)/ u s, while the time required for the precursor to reach that width is 
given by t acc (p^) oc ld(p^)/u s (see eq. [9]). Hence we choose diffusion length and time for p\ 
x = k/u and i = k/u 2 , with u = u s ^ and k = n(p ,p^), as the normalization units for length 
and time. For the normalization units for mass, we choose p = po- Then the normalized 
quantities become x = x/x, t = t/t, u = u/u, k = k/k, and p — p/p. In addition, the 
normalized pressure is expressed as P = P/(pu 2 ). With these choices, we expect that at 
time t ~ 1, the precursor width would be x ~ ldij>^) ~ 1, for example. It should be clear 
that the physical contents of our normalization are ultimately determined by the value of 
p^ anticipated to correspond to t ~ 1 as well as by the form of K,(p,p). In the simulations 
reported here, p^ was selected to give us the maximum span of p that is consistent with our 
ability to obtain converged results with available computational resources. Our choice of p^ 
is especially dependent on the nonrelativistic momentum dependence of n(p). In particular, 
when the dependence is steep, K(pmj) and ld(Pmi) can become extremely small compared to 
their relativistic values, necessitating very fine spatial resolution around the subshock. 

In Table 1, we list our numerical models classified by To and k. For example, T6P1 
model adopts T = 10 6 K and k p i with a = 1 and v = 0, while T4Bd model adopts T = 10 4 
K and Bohm diffusion, with v — 1. In the power law diffusion models of T6P1 and 
T6Pld, ~ 10 6 is chosen for the normalization, so that k(p = 1) = k*p = 10~ 6 p. For the 
Bohm diffusion models, T6Bd and T4Bd, on the other hand, p^ ~ 10 2 is chosen, because 
the steep nonrelativistic form of the diffusion makes those models too costly for us to follow 
evolution to much higher CR momenta. 

A specific example can clarify the application of these simulations to real situations. Let 
us consider a shock with u s ^ = 1.5 x 10 3 km s _1 propagating into the interstellar medium 
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with B = 5 pG. Then in the Bohm limit that the relativistic CR scattering length equals the 
gyroradius, k* = m p c 2 /(3eB) = 6.3 x 10 21 cm 2 s _1 . For the T6P1 model, for instance, the 
normalization constants are u — 1.5 x 10 3 km s" 1 and k = 6.3 x 10 27 cm 2 s _1 , so x = 4.2 x 10 19 
cm and i = 2.8 x 10 11 s. 

On the other hand, the time evolution of these shocks becomes approximately self- 
similar, as we will demonstrate. In that case the normalization choices above are entirely for 
the convenience of computation. We will eventually replace even these normalized physical 
variables with dimensionless similarity variables. To simplify the notation in the meantime, 
we hereafter drop the tilde from the normalized quantities as defined above. 

Our simulations start with a purely gasdynamic shock of Mo = 10 at rest at x = 0, 
initialized according to Rankine-Hugoniot relations with uq — — 1, po — 1 and a gas adiabatic 
index, 7 9 = 5/3. So the initial shock speed is u S j = 1 in code units. There are no pre-existing 
CRs, i.e., P c (x) = at t = 0. 



2.3. Thermal leakage and Alfven wave transport 

Although the shock Mach number is the key parameter that determines the evolution 
of CR modified shocks, the thermal leakage injection and the Alfven wave transport are 
important elements of DSA. They were discussed in detail in previous papers including Paper 
I. So here we briefly describe only the central concepts to make this paper self-contained. 

In the CRASH code, the injection of suprathermal particles via thermal leakage is em- 
ulated numerically by adopting a "transparency function", r esc (eB, v), which expresses the 
probability of downstream particles at given random velocity, v, s uccessfully swimm ing up- 



stream across the subshock through the postshock MHD waves (IKang et al.l 120021 ). whose 
amplitude is parameterized by e B . Once such particles cross into the upstream flow, they are 
subject to scattering by the upstream Alfven wave field, so participate in DSA. The condi- 
tion that non-zero probability for suprathermal downstream particles to cross the subshock 
{i.e., r esc > for p > p- in j) effectively selects the lowest momentum of the particles entering 
the CR population. The velocity v obviously must exceed the flow speed of the downstream 
plasma, u 2 - In addition, leaking particles must swim against the effective pondermotive 
force of MHD turbulence in the downstream plasma. The parameter, e# = B /B± used to 
represent this, is the ratio of the magnitude of the large-scale magnetic field aligned with 
the shock normal, Bq, to the amplitude of the postshock wave field that interacts with low 
energy particles, B±. It is more difficult for particles to swim up stream when th e wave 



turbulence is strong (eg is small), leading to smaller injection rates. iMalkov fc Volkl (119981 ) 
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argued on plasma physics grounds that it should be 0.25 < e# < 0.35. Our own CR shock 
simulations established that e# ~ 0.2 — 0.25 leads to injection fractions in the range of 
~ 10~ 4 — 10~ 3 , which are similar to the commonly adopted values in other models (e.g., 
Malkovl 1 19971 ; lAmato &: Blasil 120051 ) . In this study, we use e# = 0.2 for numerical models, 
although the choice is not critical to our conclusions. 

The CR transport in DSA is controlled by the intensity, spectrum and isotropy of the 



Bell 


1978; 


Lucek & Bell 


2000) 



there has been much emphasis on the possible amplification of the large-scale magnetic 



field via non-resonant wave-particle interactio ns within the shock precursor (e.g.. lBellll2004l ; 
Amato fc Blasil 120061 ; IVladmiriov et al.ll2006l ). Those details will not concern us here; we 
make the simplifying assumption that the Alfvenic turbulence saturates and that scattering 
isotropizes the CR distribution in the frame moving with the mean Alfven wave motion (see 
eq. [5]). Since the upstream waves are amplified by the CRs escaping upstream, the wave 
frame propagates in the upstream direct ion; i.e., u w > 0. Downstream , various processes 
should isotropize the Alfven waves (e.g., lAchterberg Blandfordlll986l ). so the wave frame 
and the bulk flow frame coincide; i.e., u w = 0. This transition in u w across the subshock 
reduces the velocity jump experienced by CRs during DSA. Since it is really the velocity 
jump rather than the density jump that sets the momentum boost, this reduces the accel- 
eration rate somewhat when the ratio of the upstream sound speed to the Alfven speed is 
finite. An additional effect that has important impact is dissipation of Alfven turbulence 
stimulated by the streaming CRs. That energy heats the inflowing plasma beyond adiabatic 
compression. The detailed physics is complicated and nonlinear, but we adopt the com- 
mon, simple assumption that the dissipation is local and that the wave growt h saturates , 



so that the dissipation rate matches the rate of wave stimulation (see eq. [I]) (jJoneslll993 



Berezhko fc Volkl 119971 ) . This energy deposition increases the sound speed of the precursor 
gas, thus reducing the Mach numbe r of the flow into the subshock, again weakening DSA to 
some degree (e.g.. lAchterbergill982l ). Thus, the CR a cceleration becomes less efficient, when 



the A lfven wave drift and heating terms are included (IBerezhko fc Volk!ll997l ; iKang &; Jones 
2006h . 



The significance of these effects can be parameterized by the ratio of the magnetic field 
to thermal energy densities, 9 = EB,o/E t h,o, in the upstream region, which scales as the 
square of the ratio of the upstream Alfven (va) and sound speeds. In Paper I, we considered 
0.1 < 9 < 1; here we set 9 = 0.1. The dependence of shock behaviors on that parameter 
are outlined in Paper I. The 9 parameter can be related to the more commonly used shock 
Alfvenic Mach number, Ma,o = Us,i/vA,o, and the initial sonic Mach number, Mo, as M4.0 = 
M^l g (lg -I)/ (29), where va,o = -Bo/V^Po- With 7 5 = 5/3 and 9 = 0.1, this translates 
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into Ma,o = 2.36M . So, for our M = 10 shocks, Ma,o ~ 24. Our initial shock speeds are 
u s>i = 150 km s _1 for T = 10 4 K and u S:i = 1500 km s" 1 for T = 10 6 K, corresponding, 
then, to va = 6.4 km s _1 and va = 64 km s _1 , respectively. For our example magnetic 
field, i?o = 5 pG, the associated upstream gas density would be p ~ 5 x 10~ 24 g cm' 3 and 
p ~ 5 x 10~ 26 g cm -3 , respectively. 



In the early evolutionary stage, as CRs are first injected and accelerated at the sub- 
shock, upstream diffusion creates a CR pressure gradient that decelerates and compresses 
the inflowing gas within a shock precursor. This leads to a gradual decrease of shock speed 
with respect to the upstream gas (Fig. 1 [a]-[b]). As the subshock consequently weakens, 
the CR injection rate decreases due to a reduced velocity jump across the subshock. The 
CR spectrum near p in j also steepens (Fig. 1 [c]-[d]). The total compression across the entire 
shock structure actually increases to about 5 in the Mach 10 shocks reported here. The 
highest momentum CRs respond to the total shock transition, which flattens the spectrum 
at higher momenta; i.e., the CR spectrum evolves the well-known concave curvature be- 
tween the lowest and the highest momenta. Each of these evolutionary features continue to 
be enhanced until preshock compression, CR injection at the subshock, and CR accelera- 
tion through the entire shock structure all reach self- consistent dynamical equilibrium states 
(Fig. 1 [e]-[f]). Once compression in the precursor reaches the level at which DSA begins 
to saturate, meaning the reduced subshock strength reduces CR injection to maintain an 
equilibrium, the shock compression (a s = pil P\ and cr t = P2/P0) as well as the gas and CR 
pressures should remain approximately constant during subsequent shock evolution. From 
that time on the structure of the precursor and the CR spectrum must evolve in tandem to 
maintain these dynamical features. 

The CR pressure is calculated from the particle distribution function by 



To see how P c evolves during the early, nonrelativistic stage, consider the idealized the 
test-particle case where the CR distribution has a power-law form, g{p) = go(p/Pinj)~ S up 
to p — p max , where < 5 = (4 — a s )/(a s — 1) < 0.5 for the shock compression ratio of 
4 > a s > 3. Then one can roughly express P c oc [(PmaxMnj) 1 " 5 - 1] oc (PmaxMnj) 1- ' 5 for 
Pinj Pm&x < 1- In a strong, unmodified shock, 1 — 5 ~ 1, and P c initially increases quickly 



3. 



Results 



3.1. Evolution toward an asymptotic state 




(8) 
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as P c oc Pmax/Pinj- We will show in §3.3, as the shock becomes modified toward the dynamical 
equilibrium state, that the CR pressure is dominated by relativistic particles and the CR 
spectrum evolves in a manner that leads to nearly constant postshock P c 2 . These features 
in the evolution of P c ^ are illustrated in Figure 1 (e) - (f). The time- asymptotic states are 
slightly different among different models, because the numerically realized CR injection rate 
depends weakly on n(p). 



The mean acceleration time for a particle to reach p n 
limit of DSA theory is given by (e.g., |Prurylll983l ) 



from p in j in the test-particle 



M - U 2 



Pmax 



Pinj 



«0 ^2_\ dp 

u u 2 J v 



(9) 



For power-law diffusion with density dependence, k p \ = K,*p a (p / p) u , the maximum momen- 
tum can be estimated by setting t = t acc as 



Pmax (t) 



a(o t - 1 ) u\ 



1/a 



1 1 1/° 
K* 



(10) 



where f c = a(a t — 1)/ [3cr t (l + crl^ u )] is a constant factor during the self-similar stage and 
u s is the shock speed in the time-asymptotic limit. As the feedback from CRs becomes 
important, the shock speed relative to far upstream flow is reduced, typically about 10-20 % 
for the shock parameters considered here (i.e., u s ~ [0.8 — 0.9]w S) j). With a = 1 and v — 1, 
for a typical value of a t ~ 5.3 for a M = 10 shock, f c w 0.13. 

In an evolving CR shock, at a given shock age of t, the power-law spectrum should 
extend roughly to p mauX (t) above which it should decrease exponentially. Then the diffusion 
length of the most energetic particles increases linearly with time as 



^max(^) 



«*PmaxW 



U, 



f C U s t. 



So l max (t) depends only on the characteristic length u s t, independent of the size of the 
diffusion coefficient, although at a given time the particles are accelerated to higher energies 
with smaller values of n*. Since the precursor scale height is proportional to / max , the 
precursor broadens linearly with time, again independent of the size of k*. This is valid even 
for the Bohm diffusion if p max 3> 1, since Kg ~ K*p for p ^> 1. Thus, the hydrodynamic 
structure of evolving CR shocks does not depend on the diffusion coefficient, even though 
the CR diffusion introduces the diffusion length and time scales in the problem. 



-12- 



3.2. Shock structure and CR spectrum in self-similar stage 

After the precursor growth reaches a time-asymptotic form, the shock structure follows 
roughly the self-similar evolution and stretches linearly with time, as noted above. Thus, 
we show in Figure 2 the evolution of a M = 10 shock with T6Pld model in terms of the 
similarity variable, £ = x/(u s ^t), for t > 1 (i.e., later stage of the shock shown in Fig. 1). 
The time-asymptotic shock speed approaches u s = uq ~ 0.9u St i for these shock parameters. 
The reduction in shock speed results from the increase in <r t , so depends upon the degree of 
shock modification. Here a t ~ 5.3, a — 1, v = 1, so equation (TTTT) give Z max ~ 0.13w s t, which 
corresponds to the precursor scale height in terms of £, = / max /(w Si i t) ~ 0.12. 

We also show the approximate self-similar evolution of the shock structure for four 
additional models with n(p,p) listed in Table 1 (Fig. 3). As discussed in §3.1, the overall 
shock structure at a given time t is roughly independent of the diffusion coefficient, except 
for some minor details in the shock profile that have developed in the early stage. Also the 
shock evolution seems to be approximately self-similar in all the models, as shown in the 
middle and right panels of Figure 3. Of course, with different values of k* and a, on the 
other hand, the highest momentum of the CR spectrum at a given time depends on k (see 
Fig. 5). 

Figure 4 (a)-(b) shows how the particle distribution at the subshock, g s (p) = f(x s ,p)p /L , 
evolves during the self-similar stage, extending to higher p max . For this model equation ffTUl) 
gives p max ~ (0.1//?*) t = 10 5 t. This estimate is quite consistent with the evolution of g s {p) 
shown in this figure. The peak value of g s (p) near p max seems to remain constant during the 
self-similar stage. This reflects the fact that P C 2 remains constant, as it must once DSA is 
saturated, and the fact that P c is dominated by relativistic CRs near p max for strong shocks. 

The injection momentum, p- m j oc yP^p/ 1 P2, becomes constant in time after the initial 
adjustment, because the postshock state is fixed in the self-similar evolution stage. Then 
the value of g s {p\n]) is fixed by g s ,th{.Pin}), the thermal distribution of the postshock gas at 
Pinj, and stays constant, too. 

Let us suppose particles with a given momentum p\ experience on average the velocity 
jump over the diffusion length £i = ld(pi)/(u s ti), Aw(£i), at time t±. At a later time t they 
will be accelerated to p = pi ■ (t/ti) 1 ' 01 and diffuse over the scale, £ = ld(p)/(u s t) = £i. 
So they experience the same velocity jump Aw(£i), as long as the velocity profile, w(£), 
remains constant during the self-similar stage. Then the spectral slopes plotted in terms of 
p/Pm&x should retain a similar shape over time. The slope of the distribution function at the 
subshock, q = — din g s /d\np+4:, and the slope of the volume integrated distribution function, 
Q = —dlnG/dliap + 4 (where G = f gdx), as a function of p/p ms , x (t) are shown in Figure 4 
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(d). Low energy particles near p ir ^ experience the subshock compression only, while highest 
momentum particles near p max feel the total shock compression. So q{p) ~ q s = 3a s / (cr s — 1) 
for p ~ p in j, while q(p) m q t = 3cr t /(a t — 1) for p ~ p max - The numerical results are roughly 
consistent with such expectations. 

Consequently, to a good approximation, g s {p) evolves with fixed amplitudes, g s {p-mi) 
and g s (pmax)) an( i with fixed spectral slopes, q s and q t at p in j and p max , respectively, while 
stretching to higher p max (^)- The volume integrated distribution function, G(p), also displays 
a similar behavior as g s (p)- In the bottom panels of Figure 4, G(p)/t and G(Z)/t are shown, 
noting that the kinetic energy passed through the shock front increases linearly with time. 

In Paper I, based on the DSA simulation results for t < 10, we suggested that the 
distribution function may become self-similar in terms of the momentum similarity variable, 
Z, defined in §1. If we define the "partial pressure function" as 



then the CR pressure is given by P c oc J °° F(Z)dZ. We suggested there that the postshock 
CR pressure stays constant because the evolution of F(Z) becomes self-similar. As can be 
seen in Figure 4 (b)-(c), the functions g s (Z) and F S (Z) at the subshock seem to change very 
slowly, giving the false impression that F S (Z) might be self-similar in terms of the variable Z. 
However, the constant shape of F(Z) cannot be compatible with the self-similar evolution of 
the precursor and shock profile. Since f s (p) oc (p/pi n j)~ 9s at Z ~ and f s (p) oc ( y p/p ma , x )~ qt 
at Z ~ 1 with constant values of p in j, q s , and q t , the shape of F(Z) should evolve accordingly 
in the self-similar stage (see Fig. 9 below). 

Figure 5 shows how the evolution of g s {p) depends on the diffusion coefficient and 
preshock temperature, while other parameters, M = 10, e# = 0.2, and 6 = 0.1, are 
fixed. The same set of models is shown as in Figure 3. The shape of g s {p) is somewhat 
different among different models, although it seems to remain similar in time for a given 
model. The causes of such differences can be understood as follows. First of all, the value 
of g s {pin}) ~ g s ,th(Pinj) depends on the value of pi n j oc (u s /c) oc M \/T^. Secondly, the 
numerically realized "effective" value of the injection momentum depends on the diffusion 
coefficient and grid spacing, leading to slightly different injection rates and shock structures. 
Thus the postshock P Cj2 and the compression ratios (i.e., the shock structure) depend weakly 
on diffusion coefficient, as shown in Figure 1 (e)-(f). The ensuing CR spectra have slightly 
different values of q s and q t as shown in Figure 5. 

The spectral slope of the CR spectrum is determined by the mean velocity jump that 
the particles experience across the shock structure. Here we examine how the precursor 
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velocity profile depends on the diffusion model. Figure 6 (a) shows the velocity structure 
U(£) = — tt(£) in the precursor (£ > 0) for five different diffusion models, where «(£) is 
defined as shown in Figure 2. We use the velocity data in the finest-level grid as well as in 
the base grid. The velocity profiles are quite similar in all the models except that the model 
with k oc p 1 / 2 shows a slightly different pattern at small scales (log£ < —5). 

Since the particles with momentum p feel on average the velocity jump over the corre- 
sponding diffusion length, we can find the velocity U(£ p ) at the distance from the shock that 
satisfies x = ld(j>) = £ P • (u s ^t). Using equation (fT0|) . we find then £ p = f c (u s /u S:i )(p/p max ) a . 
Then the particles with the same ratio of p/p max diffuse over the same similarity scale, £ p , 
and feel the same velocity jump, U(£ p ) + u w (£ p ) — U2 a cross the shock. Thus th e spectral 



slope can be estimated from the velocity profile as (e.g., iBerezhko fc Ellison! Il999l ) 



3(17 + u w ) d\n(U + U2) 

Quip) = q U {Q = r— r— + — . 13 

U + u w — il-i amp 

Figure 6 (b) shows the spectral slope, q u , which is calculated from numerical results of U + u w 
for different models. These curves compare to the q(p) curves in Figure 5. 

The numerical convergence issue should be discussed here. The base grid had a spa- 
tial resolution Axq = 2 x 10~ 3 in the code units. The small region around the subshock 
was refined with a number of levels increasing to eight, giving there a spatial resolution 
Ax$ = 7.8 x 10~ 6 . This structure was sufficient to produce dynamically converged so- 
lutions as discussed in Paper I. The diffusion length near p- m j pa 10~ 2 is, for instance, 
^dOinj) ~ K(p in j)/u s ,i pa 10" 8 in T6Pld model and l d (pmj) ~ 2 x 10~ 5 in T6P1/2 model, 
where all quantities are given in the code units. So the solution for equation dSJ) is not 
resolved for the lowest energy particles in T6Pld model, while it should be well resolved 
in T6P1/2 model. Since low energy particles cannot see the flow structure shorter than 
the minimum numerical thickness of the subshock, i.e., Axs, corresponding to the effec- 
tive diffusion length of p ~ 10 for T6Pld model, all particles below p < 10 feel the same 
subshock compression, independent of their diffusion lengths. This leads to a more or less 
constant q(p) pa q s for p < 10. The models shown in Figures 4 and 5 exhibit this trend 
except T6P1/2 model in which the diffusion of the injected particles are well resolved with 
Ax 8 // d (p inj ) = 0.4. 

The momentum integration of g(x,p), i.e., the CR pressure, is self-similar in the spatial 
similarity variable £. Moreover, the CR distribution at the subshock, g s (Z), and the volume 
integrated distribution, G(Z), both change very slowly in time, when they are expressed in 
terms of Z. So we expect that the distribution function g in the plane of (£, Z) should change 
only secularly during the self-similar stage, although, as mentioned before, g[Z) does not 
evolve self-similarly in the Z space (Fig. 7). The phase space distribution of g(£,Z) shows 
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that most of low energy particles (Z < 0.5) are confined within —0.2 < £ < 0.1, while the 
highest energy particles [Z ~ 1) diffuse over — 1 < £ < 1. Thus far away from the subshock, 
both downstream and upstream, relativistic particles dominate the CR energy spectrum. 



3.3. Analytic approximation for CR spectrum 



Based on the results of DSA simulations described in the previous subsections, we 
suggest that the CR spectrum at CR shocks with Mq > 10 in the self-similar stage can be 
approximated by the sum of two power-law functions with an exponential cutoff as follows: 

for p max > 1 >Pinj, 
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where q s > 4 and q t < 4. The specific functional form of the exponential cutoff was found 
by fitting the numerical simulation results (see. Figs. 4-5). We have shown that, after the 
precursor has developed fully, the CR pressure at the subshock approaches a time-asymptotic 
value, which leads to the self-similar evolution of the entire shock structure. Then the 
parameters, p in j, q s and q t as well as g ~ 9s,th{Pinj), become constant in time. Also, the 
value of gi seems to stay roughly constant, according the simulation results. We will show 
below gi has to be approximately constant, if P c ^ remains constant during the self-similar 
stage. Then the only time-dependent parameter in equation ffT4l) is p max (t), which can be 
estimated from equation ffTO]) . 

Now let us examine how P c2 evolves in time with the proposed form of g s {p) as p max 
increases to large values. Adopting a = 1, the contributions due to the low and high energy 
components can be calculated as 
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In Figure 8, we show the values of Pl/ go and Pujgx as a function of p max for several values 
of q s and q t and p in j = 10~ 2 . In M = 10 shocks the typical values of the compression ratios 
are a s ~ 3.1 and o t ~ 5.0, so q s ~ 4.4 and q t w 3.75. The plot shows that both Pi/go and 
Ph/ 1 g\ become constant as p max becomes ultra-relativistic, if the shock flow is modified so 
that <j s — > 3 and a t ^> 4. This explains why P Cj2 approaches an asymptotic value as p max 
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becomes large, leading to the self-similar evolution stage, after the subshock weakens to the 
subshock Mach number, M± ~ 3 — 4 and the total compression becomes greater than 4. 
Therefore g\ should stay constant, if P Cj 2 becomes constant in the self-similar stage. 

The amplitude g± can be estimated, if, for example, P Cj 2 is known from the DSA simula- 
tions; i.e., the CR pressure obtained with the proposed analytic form of g s should be equal to 
the value of P C 2 from the DSA simulations. Alternatively, as outlined in the appendix, em- 
pirical scaling relations established from simulations can connect P Ci 2 through simple physics 
to basic shock parameters. Then all the parameters necessary to construct approximations 
to the CR distribution function as given in equation (fT4"j) at arbitrary time t are known for 
the self-similar evolution stage. Since the time-asymptotic, self-similar solution of evolving 
CR shocks cannot be found (semi-) analytically either from the conservation equations or 
from the boundary conditions, we have to rely at least in part on numerical simulations to 
estimate the parameters p; n j, go, a s , a t , and P Cj 2 for given shock parameters. The analytic 
fitting forms that can approximate the DSA simulation results are described in the appendix. 

In Figures 4 and 5, we compare the analytic fitting formula in equation 0141) with the 
results of our DSA simulations. They show good agreement. These plots also demonstrate 
that g s (p max ), and therefore, gi, remains constant in the self-similar evolution stage. The 
compression ratios shown in Figure 1 are o~ s ~ 3.2 and o t ~ 5.0, so the power-law indices 
calculated with these ratios are q s = 4.36 and q% = 3.75. But the numerical value of 
q = —d In f s / dhxp near p in j is 4.2, because the diffusion of low energy particles is not resolved 
fully. The minimum value of q = —dhxf s /dlnp near p max is 3.79, slightly larger than q t , 
because of the exponential cutoff. Just to demonstrate how the proposed form of g s {p) 
fits the simulation result s, we use q s = 4.2 and q t = 3.76 instead for the curve shown in 



Figure 4. We note that iBerezhko fc Ellison! (119991 ) suggested the minimum value of q is 



<2min = 3.5 + (3.5 — 0.5cr s )/(2er t — cr s — 1). With our compression ratios, o~ s = 3.2 and a t = 5.0, 
this gives g min = 3.83, which is slightly larger than our estimate of 3.79. 

Using equations (fTOj) and (fl4"j) . we can estimate the CR spectrum g s at arbitrary time in 
the self-similar stage, as demonstrated in Figure 9 . Here the value of g\ is fixed by setting 
Pc,2 — 0.30 at t — 1 and then the same value of gi is used for the time t > 10. From the curves 
of cumulative F s (< Z), we can see that P c 2 stays almost constant with the constant value of 
gi, even though p max increases five orders of magnitude. In fact, P c ^/ {pou 2 si ) increases from 
0.30 to 0.32 as p max increases from 10 5 to 10 10 . For such a long span of time, however, g s (Z) 
or F S (Z) does not keep the same shape. At t = 10 5 , the maximum momentum corresponds 
to p max ~ 10 19 (eV/c) for protons. 

One might ask how we can justify the validity of the proposed form of g s at t ^> 1, 
while our DSA simulations have been carried up to t ~ 10 — 20. In the T6Pld model, 
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Pmax ~ 10 6 at t — 10. So, most CRs are already ultra-relativistic, and the CR spectrum 
evolves as expected (i.e., according to eq. [Ej). As long as P Cj2 stays constant, the self- 
similarity of the precursor/subshock structure would be preserved even for t ^> 1. The 
stretching of the u(x) profile in the precursor should influence the slope of the CR spectrum 
in a self-consistent way as shown in Figure 6. There is no physical reason why such feedback 
between the precursor structure and the CR spectrum cannot be extended to t 3> 1, as long 
as the assumed CR diffusion model remains valid and the most energetic particles remain 
contained within the system. In realistic shocks, however, the assumption for Bohm diffusion 
could break down due to inefficient generation of waves in the precursor. Moreover, highest 
energy particles escape from the system, when their diffusion length becomes larger that the 
physical extent of the shock. The effects of escaping particles will be explored further in the 
next section. 

We have focused here on moderately strong shock evolution with Mq > 10, since it 
is much more complicated to study nonlinear DSA at weaker shocks with Mq < 10. Non- 
relativistic CRs play a more significant role within those shocks. For instance, since P c is 
not dominated by relativistic CRs, we need to follow more accurately the diffusion of non- 
relativistic particles on scales close to the physical subshock thickness. Consequently, the 
diffusion model and the numerical grid resolution become important. The solutions also 
depend sensitively on the injection momentum, especially for shocks with Mach numbers, 
Mo < 2.5, where modifications are small, so the nearly test-particle CR spectrum is largely 
controlled by the injection momentum. Physics of thermal leakage injection, however, is not 
fully understood yet and we have only a working numerical model. Thus we defer discussion 
of semi-analytic discussion of evolving weak CR shocks to a separate paper. 



3.4. Steady State Shocks with a fixed p uh 

In realistic shocks, p max (£) may reach an upper momentum boundary, p uh , beyond which 
CRs escape upstream from the shock due to the diffusion length, Z max , approaching the 
physical size of the shocked system, or to lack of scattering waves at resonant scales of most 
energetic particles. From that time the precursor will cease to increase in scale and the self- 
similar evolution makes a transition into a stationary shock structure, or the one controlled 
by the overall dynamics of the situation. Because the shock energy is lost through particles 
escaping the system beyond p u b, the self-similar broadening of the precursor is replaced by 
a constant precursor structure in steady state. 

We have calculated additional runs for the T6Pld model in which an upper momentum 
boundary condition, i.e., g(p) = 0.0 for p > p u ^ is enforced. In these simulations once 
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Pmax(^) has reached the given value of p u b, the highest energy particles escape from the 
shock, the CR spectrum becomes steady and the precursor stops growing. Figure 10 shows 
the results of T61Pd model with p n ^ = 10 5 and without the upper momentum boundary. 
The distribution function g s (p) at the shock as well as the precursor and subshock structures 
all become steady after t > 1 in the run with p u b = 10 5 . In the other run without particle 
escape, the precursor continues to broaden and p maK (t) increases with time. However, the 
postshock states (e.g., p 2 and P C)2 ) in the two runs are quite similar and g s (p) in the steady 
state limit is almost the same as that of the run without particle escape at t « 1, except 
the exponential tail above p max . In Figure 8 we showed that P c2 stays constant as p max (t) 
increases with time, if g s (p) follows the form given in equation ffT4l) . This explains why P c2 
are very similar at different times in the two runs. Minor differences are slightly lower P c2 
and higher p 2 in the run with particle escape at p^. We note that the compression ratio 
greater than 4 results mainly from the combined effect of the precursor compression and the 
subshock jump, i.e., a t = a p ■ a s , regardless of particle escape. Energy loss due to escaping 
particles enhances the compression behind the shock only slightly in this shock, since the 
loss rate is not significant. 

In Figure 11 (a) and (b) snap shots are shown at t = 1 for the runs with p n ^ = 10 4 
and 10 5 , and at t — 10 for the run with p u b = 10 6 . For comparison, we also show the 
time-dependent solutions at t — 1 and 10 for the run without particle escape, since in the 
evolving shock p max ~ 10 5 and 10 6 at t — 1 and 10, respectively, for the T6Pld model. 
(At t = 0.1, p m ax would reach roughly to 10 4 , but by that time dynamical equilibrium has 
not been achieved and the self-similar evolution has not begun yet in the simulations.) The 
precursor structure shown in the profile of P c reflects the diffusion length of highest momenta, 
ld{Pub) oc p uh or Idipraex) oc p max {t). Here the CR pressure is plotted against £ = x/(u Sti t), 
since the results at two different times are shown together. So for example, the precursor 
width in £ is the same for the run with p u ^ = 10 5 at t = 1 (dashed line) and the run with 
p u b = 10 6 at t = 10 (long dashed line). Compared to these two runs, the run without 
particle escape at t — 1 and 10 (solid lines) have a wider precursor due to the particles in the 
exponential tail above p max (t). In Figure 11 (c) and (d) we demonstrate that the evolution 
of the shock structure is quite similar and the shock approaches similar asymptotic states 
for all the runs, almost independent of p u b or p max (t), which is consistent with Figure 8. 
The asymptotic value of P Cj2 is slightly lower and the precursor width is smaller in the 
runs with smaller p u b, as expected. Otherwise, the steady solutions with different p u h are 
approximately the same as the time-dependent solutions at the time t when p ma . x (t) equals 
to p u b. Thus the proposed form of g s {p) can be applied to steady state shocks with an upper 
momentum boundary p uh = p max as well, ignoring the exponential tail above p max - Even in 
the case where the shock structure is significantly affected by the energy loss due to escaping 
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particles, equation (fT4"j) can provide the steady state solution for g s (p), if the shock structures 
(a s , a t and postshock states) are known. 

4. Summary 

We have studied the time-dependent evolution of the CR spectrum at CR modified 
shocks in plane-parallel geometry, in which particles are accelerated to ever higher energies; 
that is, the maximum momentum p max is not prefixed. We adopted Bohm diffusion as well 
as the diffusion with the power-law momentum dependence of n(p) oc p° with 0.5 < a < 1. 
Thermal leakage injection of suprathermal particles into the CR population at the subshock 
and finite Alfven wave transport are included. Simulation parameters target nonrelativistic 
shocks with Mq > 10 in warm photoionized and hot shock-heated astrophysical environments 
with magnetic field strengths somewhat below equipartition with the thermal plasma. 

Unlike gasdynamic shocks, the time-asymptotic dynamical state of the evolving CR 
modified shocks under consideration here cannot be found analytically either from the con- 
servation equations or from the boundary conditions. So we rely on the kinetic simulations 
of diffusive shock acceleration to find the time-asymptotic state in the self-similar evolution 
stage. The general characteristics of the evolution of shock structure and particle spectrum 
can be summarized as follows: 

1) The width of the precursor, H, scales with the diffusion length of the most energetic 
particles and for diffusion that scales as k = n*{po/ p) v p a , increases linearly with time, 
i.e., H ps Z max ~ 0.1u s t, independent of the magnitude (ft*)and the value of a. 

2) If the acceleration time scale to reach relativistic energies from injection is much 
shorter than the dynamical time scale of the shock system (i.e., k* <C 0.1u s R, where R is 
the characteristic size of the shock), the CR pressure at the subshock approaches a constant 
value as the P c at the shock becomes a significant fraction of the momentum flux through 
the shock, ~ Po M o- For typical nonrelativistic shocks associated with cosmic structure this 
transition roughly corresponds to a time when p max becomes ultra-relativistic. Once this 
dynamical equilibrium develops, the shock precursor compression and the subshock jump 
are steady, leading to a self-similar stretching of the precursor with time. Consequently, the 
subshock compression ratio, a s , the total compression ratio, a t , as well as the postshock gas 
and CR pressures, P g ^ and P c ^, remain constant during the self-similar stage of the shock. 

3) The lowest energy particles diffuse on a = n(pmi)/u s and, so, experience 
only the compression across the subshock. Thus, near the injection momentum, p in j, the 
CR distribution function is given by f(p) ps f s ,th(Pmj)(p/Pm})~ 9s where f Sjt h is the thermal 
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Maxwellian distribution of the postshock gas and q s = 3cx s /(a s — 1). The amplitude fthipioi) 
is determined by the thermal leakage injection physics, since that establishes 



4) The most energetic particles diffuse on a scale Z max = K{p max )/u s and, so, experi- 
ence the total compression across the entire shock structure. Consequently, near p max , f(p) 
flattens to (p/p m &x)~ qt , where q t = 3a t /(cr t — 1). For p > p max , f(p) is suppressed by an 
exponential cutoff. 

Considering these facts, we proposed that the CR spectrum at the subshock for arbitrary 
time t after self-similar evolution begins can be described approximately by the following 
simple analytic formula: 
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where / = f s ,th(Pmy> and p max oc (u^t/n*) 1 ^ is given in equation (ITU1) . The parameters, 
Pinj, q s and qt can be estimated from the shock structure in the self-similar stage using DSA 
simulations results as outlined in the appendix. The amplitude, f\, has to satisfy the relation 
^(Pmax) = fs(Pmax)Pmax ~ constant in order for the postshock P c to remain steady. So, the 
momentum distribution function g(p) is shifted to higher J9 max in time, while keeping the 
amplitude at p max constant in the self-similar stage. Hence p max is the only time-dependent 
parameter in equation (TIBl) . 



In a realistic shock geometry, however, CRs may escape upstream from the shock due to 
largest diffusion length approaching the physical size of the shocked system, or due to lack of 
scattering waves at resonant scales of most energetic particles. Once p max approaches some 
upper momentum boundary at p np , the shock structure and the CR spectrum develop steady 
states that are approximately the same as the evolving forms with p max = p up , except that 
some differences in the shock structure due to energy loss from escaping particles. Otherwise, 
the shock structure parameters and the approximate analytic form for the CR spectrum in 
the self-similar stage are c onsistent with previously proposed analytic a nd semi-analytic 
steady state solutions (e.g.. iBerezhko fc Ellison! Il999l ; lAmato fc Blasil 120051 ) . 



Finally, we note that the evolution of the CR spectrum is secular in terms of the variable, 
Z = ln(p/pi n j)/ ln(p max /p in j), which alluded wrongfully the self-similar evolution of the partial 
pressure function F S (Z) in Paper I. In fact there is no similarity relation between p and t. 
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A. Analytic approximations for dynamical states 

As we noted in the Introduction, there are several analytic and semi-analytic treatments 
of strong, steady-state CR modified shocks. The full time-asymptotic state of evolving CR 
modified shocks can be obtained only through numerical simulations of nonlinear DSA. 
However, such simulations show strong similarities between steady-state and asymptotic, 
evolving shocks. Here we outline some of those basic dynamical relations as they can be 
estimated analytically and empirically from our simulations, as reported in this paper and 
previously in Paper I. 

A key to this comparison is the fact that the time scale for evolution of the shock 
precursor is the acceleration time scale to reach p max , t acc ~ 10(Z max AO (see eq. 0), which 
is characteristically an order of magnitude greater than the time scale for a fluid element to 
pass through the precursor, t dyn ~ l m &x/ u s- Then, in following a fluid element through the 
precursor, one can neglect terms d/dt compared to terms ud/dx in evaluating the Lagrangian 
time variation, d/dt. For example, equation which can be expressed as 



d f Pg\ 2 W 



dt\p 5 / 3 J 3p 5 / 3 ' 
assuming 7 9 = 5/3, then gives for an evolving precursor 
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where <r p = p±/po is the precursor compression factor. The quantity 
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dx (A3) 



was introduced in Paper I, and measures entropy added by Alfven wave dissipation while the 
fluid element crosses the precursor, normalized by MqPo/Po • Since equation (1A2j) applies 
to an evolving shock, the subscripts '0' and '1' refer to states of a given fluid element as 
it enters the precursor and as it reaches the subshock. The approximation comes from 
neglecting explicit time variations in \W\ and p in evaluating /. Equation (1A2j) is exact for a 
steady state shock. In the absence of Alfven wave dissipation, this equation simply states the 
properties of adiabatic compression through the precursor, which obviously does not depend 
on the precursor being steady state. 
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Along similar lines, momentum conservation of a fluid element passing through the 
(slowly) evolving precursor gives 



PcA + PqA ~ -Pg,0 



1 

f 



(A4) 



which can be combined with equation (IA2I) to produce a simple estimate for the CR pressure 
at the subshock, 



PcA = Pc 



c,2 



Poul 



1 



<7„ 



o 5/3 1 
3 0V — 1 



M 2 



5 p 



(A5) 



By substituting equation (IA5I) into equation flA3j) along with equation (4), one can obtain 

(A6) 



, 5 v A , P cl 
I « -- 



3 m Po«o' 

where, once again, the approximation reflects neglect of explicit time variation in the shock 
structure during passage of a fluid element through the shock. Substituting this back into 
equation (IA5B we obtain 



-r C ,2 



o 5/3 1 
3 ay — 1 



0V, 



M 2 



2^0 5/3 

3 n p 



(A7) 



Given P Cj i = P c>2 from equation (1A7jl and using equation (1A4I) it is straightforward to 
determine, as well, P 9t i. 

Although we can estimate approximately the postshock pressures, P g ^ and P c ^, for a 
given value of precursor compression, we must rely on numerical simulations to obtain the 
value of a p for different model parameters. In the remainder of this appendix we present some 
practical expressions for the shock dynamical properties obtained in our DSA simulations 
using a wide range of Mach numbers for the thermal injection parameter e# = 0.2, the Alfven 
wave transport parameter, 9 — 0.1 and the diffusion coefficient, k = K*p(p/p ). In Figure 11 
the time-asymptotic values of postshock CR pressure, gas pressure and compression ratios 
are plotted against the initial shock Mach number (M > 1.5). 

For M < 2.5, the CR modification is negligible, so the postshock gas pressure and the 
shock compression ratios a t = a s are given by the usual Rankine-Hondo relation for pure 
gasdynamic shocks. 

For Mq > 2.5, the numerical results for the postshock gas pressure can be fitted by 

-0.4 



P 



S,2 



PoKi 



0.4 



(Mo 



(A* 
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The time-asymptotic density compression ratios can be approximated as follows: 

f M \ ' 17 

a s » 3.2 f j for 2.5 < M < 10, (A9) 
M X 04 



a s « 3.2 ( ) for M > 10, 



a * ~ 5 '° v iT J for 2 ' 5 - M ° - 10 ' ( A1 °) 
^ ^ 5 ( T(T J forM ° >10 - 

We note that the subshock compression depends only weakly on M , while the total com- 

1/3 

pression increases approximately as M . Even for strong shocks with Mo up to 100, the 
total compression ratio is less than 10, because the propagation and dissipation of Alfven 
waves upstream reduces the CR acceleration and the precursor compression. 

The postshock CR pressure can be fit empirically as follows: 

-^§- « 2.34 x 10~ 2 (M - l) 3 for 1.5 < M < 2.5, 



Pou s>l 

\ 2 0.58(M -1) 4 2.14(M -1) 3 13.7(M -1) 2 



poul Ml Ml M 4 

27.0(M - 1) 15.0 



(All) 



M 4 M 4 



+ for 2.5 < M < 100, 

Pel 



Pouli 



0.55 for Mo > 100. 



These fits are plotted in solid lines in Figure 11. Since <j v = a t /a s , equations (A9) and (A10) 
can be used along with equation (1A7|) to estimate P C]2 (dotted line in Fig. 11). 



In iKang et al.l (120021 ) we showed that the effective injection momentum is Pi n j/pth ~ 2.5 
for Mo > 10 for the injection parameter e# = 0.2, where pth = 2 y/ kT 2 / m p c 2 and T 2 = 
(Pgfil P2){jn P /k) is the postshock gas temperature. Then the thermal distribution at the 
injection momentum, g.s,th{Pmj), can be calculated from the Maxwell distribution, since the 
postshock gas states, T 2 and p 2 , are known. 
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Table I. Preshock Temperature and Diffusion Coefficient in Numerical Models 



Model Name 


T 




k/k 


Description for Diffusion Coefficient 


T6Pld 


10 6 K 




10~ 6 p(p /p) 


power-law diffusion with p' 1 dependence 


T6P3/4d 


10 6 K 




io-V /4 (po/p) 


power-law diffusion with p~ x dependence 


T6Pf 


10 6 K 






power-law diffusion 


T6Pf/2 


10 6 K 




1.78 x 10- V /2 


power- law diffusion 


T4Pld 


10 4 K 




10" 5 p( Po /p) 


power-law diffusion with dependence 


T6Bd 


10 6 K 


10" 


-V/Vp 2 + l(Po/p) 


Bohm diffusion with p~ x dependence 


T4Bd 


10 4 K 


10" 


-V/Vp 2 + l(Po/p) 


Bohm diffusion with dependence 
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Fig. 1. — (a)-(d): Snap shots of M = 10 shock of T6Pld model up to t = 1 in terms 
of a similarity variable £ = x/(u s ^t). The flow velocity, CR pressure, and CR distribution 
function at the subshock, g s {p) = / s (p)p 4 , and g s (Z), are shown at t — 0.1 (dotted lines), 
0.2 (dashed), 0.3 (dot-dashed), 0.6 (dot-long dashed), and 1.0 (solid). The long dashed lines 
show the initial shock structure, (e)-(f): Time evolution of the postshock CR pressure and 
the compression ratios for M = 10 shocks with four different models of diffusion coefficient 
n(p) (see Table 1). 
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Fig. 2. — Self-similar evolution of Mq = 10 shock of T6Pld model. The shock structure is 
shown at t = 2 (dotted lines), 10 (dashed), and 20 (solid) as a function of the similarity 
variable £ = x/(u s ^t) in the shock rest frame. The long dashed lines show the initial shock 
structure. 
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Fig. 3. — Self-similar evolution of Mq = 10 shock of four different models listed in Table 1, 
shown at t — 1 (dotted lines), 5 (dashed), and 10 (solid). The CR pressure is shown as a 
function of x (left panels) and the spatial similarity variable £ = x/{u s ^t) (middle panels). 
The gas density is shown at the right panels. 



-30- 



o 



-6 
-7 
8 
-9 



- 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 

M =10, T6Pld 


, i 1 , 


ii- 






< *^~\ 
• \\ 






t= 2, 10, 20 


: y 

\ \ 






1 , , , 1 , , , 1 , , , 1 


': 1 
, , ; l ,1 





■2 2 4 6 8 
log(p/m p c) 



— i -7 
° -8 



~i — i — i — i — I — i — i — i — i — I — i — r 



p. .= 10" a , p =10 5 * 

rinj ' rmax 



t ^ 




q s = 4.20, q t =3.76 

I i i i i I i i i i I 







0.5 



1 



3 



O 




-1 

-2 



^ 1 II II 1 

r f. . 


I i M i ii in 

( 


i i _ 
c)j 




F S (<Z) 




1- \A i 1 i i i 1 i 


i i 1 i i i 1 i i i 1 


ti i - 



0.2 0.4 0.6 0.8 1 1.2 
Z 




-8 -6 -4 -2 
log(p/p max ) 





-6.5 


_ i | i i i 1 i 


i | i i i | i i i | i 


1 1 1 - 




-7 






(e)j 






A 






-7.5 


















-8 




















-8.5 










o 


- 1 1 i i i 1 i 


1 1 1 1 1 1 1 1 1 1 n 


, , \- 





-6.5 
-7 


: 1 1 1 


i i | i i i i | 


i i _ 

(0 j 


[G/t 


-7.5 










-8 










-8.5 

-9 


'- 1 , , 


i i 1 i i i i 1 


1 1 1 " 



-2 2 4 6 8 
log(p/m p c) 







0.5 



Fig. 4. — CR distribution function for Mq = 10 shock of T6Pld model, shown in Fig. 2, at 
t = 2 (dotted lines), 10 (dashed), and 20 (solid), (a)-(b): The distribution function at the 
subshock, g s (p) and g s (Z), where Z = ln(p/p in j)/ ln(p max /p in j). (c): The partial pressure, 
F s , defined in equation (fT2l and its cumulative distribution, F s (< z). (d): The power-law 
slopes, q = — din g s /d\np + 4 and Q = — d\nG/d\np + 4. (e)-(f): The volume integrated 
distribution function G = J gdx plotted against log(p) or Z. The long dashed lines in (a), 
(b), and (d) show the analytic fitting given in equation ( |T4|) with q s = 4.20, q t = 3.76, 
Pinj = 10~ 2 , and p max = 10 5 t at t = 10. 
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Fig. 5. — CR distribution function at the subshock, g s (p) and g s {Z), and the power-law 
slopes, q = —d\Tig s /d\iap + 4: are shown at t — 1 (dotted lines), 5 (dashed), and 10 (solid) for 
the four diffusion models shown in Fig. 3. The long dashed lines show the analytic fitting 
given in equation (fl4"l) at t = 10. The adopted values of q s and q t are given for each model. 
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Fig. 6. — (a) Velocity profiles in the precursor as a function the similarity distance from 
the subshock for five different models with M = 10 listed in Table 1. (b) The power slope 
calculated with equation (TIB"]) using the velocity profile shown in (a). 
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Fig. 7. — Left panel: Contour plots of g(£,Z) at t = 10 (dotted lines), 20 (solid lines) for 
T6Pld model. Right panel: Contour plots of g(£, Z) at t — 5 (dotted lines), 10 (solid lines) 
for T6P1/2 model. 
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Fig. 8. — Left panel: Pi, defined in equation ( fT5l) . for q s = 4.0 — 4.6. Right panel: Ph, 
defined in equation (fT5l) . for q t = 3.0 — 4.0. Here the injection momentum is p m j = 10 -2 . 
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Fig. 9. — Upper panels: CR distribution at the subshock calculated using the analytic fitting 
formula in equation (jHj) with q s = 4.20, q t = 3.76, p in j = 10~ 2 , and p max = 10 5 t. Lower 
panels: Power-law slope of the fitted g s , i.e., q = —dlng s /dlnp + 4, and partial pressure, 
F S (Z), and its cumulative distribution, F s (< Z). 
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Fig. 10. — Comparison of the run with particle escape at p u b = 10 5 (solid lines) and the 
run without particle escape (dashed line) for T6Pld model. The shock structure and the 
CR distribution are shown at t — 1 and 4. 



-37- 



0.3 



0.2 



0.1 







-III 
_ (a) 
-t=l 


i i i 1 i i _ 












-6.5 


-t=l 










4 = 10 








-7 


-t=l, 10 






m 

M 






o 






K ~- 




-7.5 


- i i i 1 i i i 1 






-8 



0.4 -0.2 
x/(u s>i t) 



i i i I i i i I i i i I i i i I i i 
(b) 



P ub =10 4 . t=l 
" p ub =10 5 , t=l 
p ub =10 6 , t=10 

-l p - ( l l) ', ^ =1 ',I 1Q 

I I I I I I I I I I I I I I I I 



2 2 4 6 
log(p/mc) 



0.3 



0J w 

0.2 



cv 
d 



0.1 







i i i mill — i i i inii — i i i Mini — i i i linn — r 



(c) 




0.001 0.01 0.1 1 10 

t 



i i i hum — i i 1 1 him — i i i Mini — i i i mill — r 




0.001 0.01 0.1 1 10 

t 



Fig. 11.— Comparison of the runs with and without particle escape at p u b for T6Pld 
model, (a) CR pressure profiles in the three runs with p u ^ = 10 4 at t = 1 (dotted line), with 
p u b = 10 5 at t — 1 (dashed), and with p ub = 10 6 at t — 10 (long dashed line). The solid lines 
are for the run without particle escape at t — 1 and 10. (b) CR spectrum at the subshock. 
(c)-(d): Time evolution of the postshock CR pressure and the compression ratios. The same 
line types are used in all the panels. 



-38- 




Fig. 12. — Time-asymptotic values of postshock gas and CR pressures in units of initial 
shock ram pressure (left panel), subshock compression ratio (triangles, right panel) and 
total compression ratio (circles, right panel) as a function of initial shock Mach number Mq 
for T6Pld models. The solid lines show our fitting formulas given in equations (1A8I) -( IA~TTT) . 
The dotted line shows the estimate given in equation (1A7I) . adopting the numerical values of 
= °~t/°~s in equations flAlOl) and flAllj) . 



